Phase diagram of a polydisperse soft-spheres model for hquids and colloids 
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The phase diagram of soft spheres with size dispersion has been studied by means of an optimized 
Monte Carlo algorithm which allows to equilibrate below the kinetic glass transition for all sizes 
distribution. The system ubiquitously undergoes a first order freezing transition. While for small 
size dispersion the frozen phase has a crystalline structure, large density inhomogeneities appear in 
the highly disperse systems. Studying the interplay between the equilibrium phase diagram and the 
kinetic glass transition, we argue that the experimentally found terminal polydispersity of colloids 
is a purely kinetic phenomenon. 



PACS numbers: 64.60.Fr,64.60.My,66.20.-|-d 

The equilibrium phase diagram of dense classical fluids 
is far from being fully understood, especially as regards 
the influence over the freezing transition of some disorder 
in the parameters of the interaction. While fluids made of 
identical atoms crystallize easily upon lowering the tem- 
perature or increasing the density, the fate oi polydisperse 
systems (colloids, for instance), where the particle size a 
is sampled from a probability distribution, P{cr), is still 
a matter of debate. 

On the theoretical side, the effect of size dispersion, 
measured by the quantity d (the ratio among the stan- 
dard deviation and the mean of -P(cr)) over the phase 
diagram has been addressed mostly in the case of the 
hard-spheres model, which is meant to describe colloidal 
systems [Ij. At least for small polydispersity, d seems 
to be the only feature of P{a) that controls the physical 
results. Different theories predict conflicting scenarios in 
the region of large d. The moment free-energy method pi] 
suggests phase separation between many different crys- 
tal phases, each one with a much narrower size dispersion 
than S (a phenomenon termed fractionation). However, 
a density functional analysis Q predicts the existence of 
a terminal polydispersity St beyond which the formation 
of the crystal is thermodynamically suppressed. Numeric 
simulations of the hard sphere model find some agree- 
ment with both the fractionation 0, and the terminal 
polydispersity [6] scenarios. As regards models with a 
soft potential (e.g. Lennard- Jones), which are more ap- 
propriate to describe liquids, no extensive study on the 
effects of polydispersity has been performed so far. 

On the experimental side, crystallization of very vis- 
cous colloidal samples with S > St ~ 0.12 does not occur, 
even after months spent from the preparation Q . Further 
evidence supporting the terminal polydispersity scenario 
comes from the finding of reentrant melting (crystal to 
fluid transition when increasing the density) on polydis- 
perse colloids in confined geometry [§]. 

Yet, these experimental results do not necessarily re- 
veal features of the phase diagram. Indeed, the processes 



needed to keep the system in thermodynamic equilib- 
rium often become exceedingly slow. Such processes are 
the nucleation of the solid within the fluid and the a- 
relaxation in the super-cooled fluid 0]. 

Here we obtain the equilibrium phase diagram of poly- 
disperse soft-spheres in the (N, V, T) ensemble, aiming 
to describe both liquids and colloids. This is made pos- 
sible by the combination of an optimized Monte Carlo 
(MC) method (which, unlike standard MC, equilibrates 
below the kinetic glass temperature) and a novel Finite- 
Size Scaling study of the fluid-solid coexistence line. At 
all 5, a flrst-order freezing transition separates the fluid 
phase from the low temperature solid. This rules out 
the thermodynamic origin of the terminal polydispersity 
scenario. However, we show that a Brownian dynamics 
will not find crystallization for 5 > 0.12, in quantita- 
tive agreement with colloids experiments 7]. Further- 
more, depending on polydispersity the solid phase can 
be either a crystal or an inhomogeneous phase (hereafter 
Tphase). The freezing temperature shows a reentrant 
behavior when increasing d. Interestingly, in the range 
6 e [0.12, 0.38] the kinetic glass transition occurs for tem- 
peratures, T, and densities, p, in the stable fiuid phase. 

Specifically, in our simulations soft-spheres of radius at 
and (Tj at distance r interact via the pair potential [sol ] 

The effect of T and p is encoded in the single thermo- 
dynamic parameter F = pT~-^/^. Although ^ general- 
izes well known models for simple liquids [IQ] , its scale- 
invariant form suggests that it could describe as well col- 
loids, whose size is in the micrometer range. We consider 
a flat size distribution, constant in the range [tTmin, Cmax]- 
In order to eliminate sample-to-sample fluctuations we 
follow fllj, which for a flat distribution amounts to pick 

(Ti = CTmin + '^{i - 1), with A = ((Tmax - 0-,„in)/ (A^ - 1). 
Note that CTmax/Cmin ^ oo at l5oo = l/\/3- 

The numerical investigation of equilibrium properties 
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FIG. 1: The F — S phase diagram shows three phases: a fluid 
at high temperatures, a crystal, and a inhomogeneous soUd 
(I-phase). The boundary between the crystal and the I-phase 
Ues at 0.18<(5c<0.21. The vertical line in the fluid phase 
is the kinetic glass transition. (Inset) Integrated relax;ation 
time [l^ . r, for e (full) and J- (open), both for standard 
(squares) and local swap (circles) MC, for A'^ = 256. The 
kinetic glass transition arises when r ~ 10® standard MC steps. 

of fluids is limited by the practical impossibility to equi- 
librate when either the relaxation time r within the 
fluid phase, or the freezing time th SlJ become com- 
parable with the time scale of the simulation. In or- 
der to achieve equilibrium we straightforwardly adapt to 
model ([1]) the local swap MC algorithm 1^, [l3| , which 
outperforms 13l . Il4l other algorithms, such as the non- 
local swap MC [iSj and the density of states MC 
In fact, the standard method of studying a first-order 
transition in a system of finite size, L, looks for a double 
peak in the histogram of the internal energy per particle 
e, accompanied by a peak of the specific heat Cy 17. 1^. 
This procedure is extremely demanding on the quality of 
the statistical sampling, and has been scarcely used (if 
at all) for glass-forming liquids or colloids. Fortunately, 
the local swap MC allows us to employ it. 

The thermalization issue needs to be addressed in three 
different regimes, see Fig. [1] the liquid phase, at the 
phase coexistence line, and the solid phase. The swap 
algorithm avoids the cage effect, thus equilibrating eas- 
ily the whole liquid phase (it reduces r by two orders of 
magnitude as compared with standard MC, see Fig. [1]- 
inset). Thermalization in the deep solid phase, not at- 
tempted in this work, would require a different approach. 
At the phase coexistence line, our criterion for thermal- 
ization required the observation of tenths of back and 
forth tunneling events between the liquid and solid phase. 
The difficulty for meeting the criterion grows dramat- 
ically with the number of particles (^ exp[i7iV^/'^], S 
being the liquid-solid surface tension). A stronger, more 
quantitative check is the consistency of the histogram 
reweighting fs^- We equilibrated N = 256 particles for 
i5 > 0.1 (and for the monodisperse system) and iV = 500 
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FIG. 2: Equilibrium quantities for 5 — 0.24 {top-left) Spe- 
cific heat vs. F, from two simulations, one at the size depen- 
dent critical point (full symbols denote the actual simulations, 
full lines coming from histogram reweighting), the other at a 
higher F (open symbols, and dashed lines). The maximum of 
the specific heat scales with A^. [Top-right) {^)/N vs. F for 
the same simulations of top-left panel. (Bottom-left) Normal- 
ized histogram, /, of the energy of the inherent structures, for 
N = 256 (empty) and N = 500 (filled), at the F value where 
the peaks of the instantaneous energy have equal height fl8|. 
Inherent structures were found every 10^ MC steps. (Bottom- 
right) As Bottom-left, for JT. In the fiuid phase is 0(1) (in 
the solid, JT is 0(N)). We also show (dashed line) eia and 
T histograms for A'^ = 864, where only two back and forth 
tunneling events between the liquid and solid phase were ob- 
served. Yet, data nicely agree with the predicted A'^-scaling. 



for 6 > 0.21. However, for local swap MC tfr remains be- 
tween 10^ and 10^ MC steps for all 6, while for standard 
MC it grows beyond 10^ steps for 6 > 0.12. 

To detect the possible existence of large density fluctu- 
ations we focus on S(q) = j exp[iq • (r^ — rj)] (ri is 
the position of the i-th particle). Note that {S{q)) is the 
static structure factor. The longest wavelengths fluctua- 
tions are studied through T = [5(^, 0, 0) -)-5(0, ^,0) -I- 
5(0, 0, ^)]/3 . In fact, in an homogeneous liquid or crys- 
tal phase we expect T to be of order 1, but for a macro- 
scopically inhomogeneous phase it becomes of order N . 

In Fig. [2] we show as an example the results obtained 
for 5 = 0.24 where an evidence for a freezing transition 
is presented. The specific-heat displays a peak of height 
proportional to N while, as usual for small systems, its 
position suffers a strong finite size shift. We found con- 
venient to use the histogram of the energy per particle of 
the inherent structures ejg, i.e. the minima of the poten- 
tial energy [19|]. The advantages are twofold (Fig. [2]): it 
largely absorbs the effects of the finite-system shift of the 
critical temperature (so that the position of the peaks is 
almost TV independent) and it makes each peak narrower. 

Note that (5 = 0.24 is much higher than the terminal 
polydispersity i5t reported in experiments and simula- 
tions [71, |n| . Actually, a freezing transition arises in the 



full range of S studied. The line of phase-coexistence, 
as located by the arising of a double peak for N = 256, 
between the fluid and the solid phase is shown in Fig. [T] 

For S = 0, we find a body-centered cubic (bcc) crystal, 
as expected for modest N [23|. Indeed, S{q) displays 
peaks of order N at q ^ ^(1, 1, 0), ^(0, 1, 1), ^(1, 0, 1) 
with lattice spacing a ^ 0.78. The same Bragg peaks 
are found at 5 = 0.12, broadened due to disorder j33[ (as 
compared with (5 = 0, the intensity reduces by a factor 
1/4, but it still increases linearly with TV). 

Interestingly, for 6 = 0.24 the histogram of J-" (Fig. ^ 
develops a double peak at the freezing point. Left-peak's 
position is iV- independent, as expected for the liquid 
phase (see also the scatter plot in Fig. [3]) . On the other 
hand, the second peak (corresponding to the solid) shifts 
right proportionally to N, revealing that the solid phase 
is not a standard crystalline one. For A^ = 500 and 6 < 0.2 
it is risky to use histogram techniques, as the tunneling 
from solid to liquid is harder to observe than the reverse 
liquid to solid tunneling, thus we show in Fig. [3] a scatter 
plot. It is clear that in the crystals we found at (5 = 0.18 
(corresponding to the lower values of eis in Fig. [3]) 
takes a A^-independent value 3J|. Summarizing, at high 
polydispersities we have found a freezing transition from 
the liquid to a solid inhomogeneous phase, while at low 
polydispersities the low temperature high density phase 
(large F's) is a standard homogeneous crystal. The study 
of the transition between the crystal and the I-phase is 
left for future work. To gain some insight about the I- 
phase we measured the intensity of the Bragg peaks for 
the inherent structures in the solid phase. More precisely, 
we define B as the maximum 35] of S{q). In a crystal 
B is of order N (the corresponding q is in the reciprocal 
lattice), while in a fluid B is always of order 1. In Fig. [4] 
we show the histogram of B along the phase coexistence 
line both for A^ = 256 and A^ = 500. For every S we find a 
double peak structure. The solid's peak shifts right pro- 
portionally N. Thus (B) and (JF) provide a classification 
of the r — 5 plane: (B) distinguishes the solid from the 
fluid phase, while (JF) is of order N only in the I-phase. 

The solid peak in Fig. |4] shifts left with increasing S 
(at S = 0.45 the two peaks are hardly resolved). The 
I-phase might signal either fractionation [1, 0| or phase 
coexistence of a crystal with an amorphous solid (Fig. 21 
top-left). Further work is needed to clarify this point. 

The existence of a freezing transition for all S rules out 
the terminal polydispersity scenario. But even if a sta- 
ble solid phase exists for large S, it might be dynamically 
inaccessible on experimental time-scales. In order to de- 
tect the occurrence of the kinetic glass transition on the 
S ~ r plane we adopted standard MC simulations which 
(at variance with local swap MC), mimic the real dynam- 
ics of the system by modeling it as a Brownian motion. 
Note that the correspondence between the single step of 
standard MC and the physical time units widely differ for 
liquids and colloids, being ^ 10~^^ sees, in the former 
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FIG. 3; Scatter plot of !F vs. eis below (bottom) and above 
(top) the I-phase-crystal transition line at the liquid-solid 
phase coexistence, for A*' = 256 (left) and A'^ = 500 (right). The 
number of points is ~ 45000 {N = 256) or ~ 15000 (Af = 500). 
The darker the shade of gray, the higher the density of points. 
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FIG. 4: (Top-left) A 5 = 0.3 sohd configuration of 500 
particles, at phase coexistence, as projected into the XY 
plane (circles are centered at particles positions, their diam- 
eter being proportional to particle sizes). Larger particles 
form the crystalline central band (the Bragg peaks's analysis 
suggests a FCC ordering). Normalized histograms of B for 
5 = 0.24,0.3,0.45 with 256 (empty) and 500 (filled) particles. 



case and ~ 10^^ sees, in the latter one 21j . 

We have located the kinetic glass transition in Fig. [1] as 
the point where r reaches 10^ (standard) MC steps. In 
the case of colloids, this roughly corresponds to a relax- 
ation time of three hours, hence our kinetic glass transi- 
tion corresponds to the experimental one. On the other 
hand, since in the case of liquids 10^ MC steps ^ 10~^ 
sees., our kinetic glass transition is rather close to the 
mode-coupling transition [22|. Indeed, for most molecu- 
lar and polymeric glass-formers t at the mode-coupling 
transition lies between 10~^'^ and 10~^'^ sees. [23*1. By 
the way, since r grows by many order of magnitudes in 
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a narrow range, for liquids the difference in F between 
the experimental and the mode coupling transition will 



not be larger than 5% [2J|. In the inset of Fig. [T] we 
observe that the kinetic glass transition is at rj, ~ 1.45 
for (5 = 0.24. This value coincides with the one found for 
the soft-spheres binary mixtures (i.e. a bimodal -P(o")) 
with (5 = 0.09 [H and (5 = 0.16 [3]. This suggests that 
Fc is independent on 6. The vertical line at Fc — 1.45 
intersects the freezing transition line at two intersection 
points, creating a sort of no man's land where equilibrium 
is experimentally unreachable on human time-scales. In 
fact, the freezing time depends on 6, growing tremen- 
dously along the phase coexistence line. When approach- 
ing the first intersection point {S « 0.12, see Fig. [T]), for 
standard MC we only know that ift. becomes larger than 
10^ steps. For 5 > 0.12 our standard MC simulations 
did not find the solid phase. We thus argue that the 
tremendous slowing down both of r and ift. makes the 
solid phase for (5 > 0.12 experimentally unattainable. 

We studied the equilibrium phase diagram in the S — F 
plane for a polydisperse soft sphere model with empha- 
sis on the first order freezing transition. Optimized MC 
algorithms and Finite Size Scaling analysis were crucial 
to obtain equilibrium properties. We found two different 
solid phases, a standard crystal at small S and a highly 
inhomogeneous phase at large 6. Moreover, for 6 lying be- 
tween 0.12 and 0.38 the glass transition is no longer pre- 
empted by the freezing transition, making this range very 
promising for the theoretical study of glasses. In experi- 
ments non-equilibrium effects should dominate the whole 
region S > 0.12, where either the kinetic glass transitions 
or the growth of the freezing time-scale are expected to 
hide the freezing transition. This is in quantitative agree- 
ment with the findings for colloids TJ suggesting that the 
polydisperse soft-spheres model ^ could catch the fea- 
tures both of molecular glass-formers and of colloids. 
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